#### "Gender Inequality and Authoritarian Regimes" ####
# authors: "Lars Pelke"
# date: 2020-01-20

## Set WD ##

## use 01_data_wrangling before running this file  OR

## load data ##

vdem <- readRDS("calculations/data/vdem_wth_merged.rds")


#### Figures A1 Main Paper ####

vdem$wth_regime <- as.factor(vdem$wth_regime)
vdem <- vdem %>%
  drop_na(wth_regime)

ggplot(data = vdem , aes(x = wth_regime, y =v2xpe_exlgender)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Political Exclusion by Gender")

ggsave("calculations/Output/robustness_wth/Figure1.pdf", height = 14, width = 24, units= c("cm"))


#### Figure 2 Main Paper ####

ggplot(data = vdem , aes(x = wth_regime, y =v2xps_party)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Party Institutionalization")

ggsave("calculations/Output/robustness_wth/Figure2.pdf", height = 14, width = 24, units= c("cm"))

#### Prepare Data for analysis ####
library(scales)

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(1,0))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(1,0))
summary(vdem$gender_inclusion_output_mean)

vdem <- vdem %>%
  drop_na(cown, year, wth_regime)

vdem2 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)


#### Descriptive Tables ####

#functions 
mean.new <- function(v) {
  if (all(is.na(v))) { return(NA) } else { return(mean(v, na.rm=T)) }
}

## Table C1 ##
table1 <- vdem %>%
  group_by(wth_regime) %>%
  summarize(mean_input = mean.new(gender_inclusion_input_mean), 
            mean_output = mean.new(gender_inclusion_output_mean))%>%
  rename("Average WPI on\nthe input side" = mean_input, 
         "Average WPI on\nthe output side" = mean_output, 
         "Regime Type" = wth_regime)

library(stargazer)
stargazer(table1, summary=FALSE)


## Table C2 ##

summary(vdem$wth_regime)

vdem <- vdem %>%
  ungroup() %>%
  mutate(wth_regime_num = case_when(wth_regime == "military" ~ 1, 
                                    wth_regime == "monarchy" ~ 2,
                                    wth_regime == "multi-party"~ 3, 
                                    wth_regime == "no-party"~ 4, 
                                    wth_regime == "one-party"~ 5)) 

vdem <- vdem %>%
  group_by(country_name) %>%
  mutate(group_id = cumsum(wth_regime_num != lag(wth_regime_num, default = FALSE))) # group_id for each regime spell

vdem <- vdem %>%
  mutate(regime_type_id = str_c(country_name, group_id, sep = "_"))

table3 <- vdem %>%
  group_by(regime_type_id) %>%
  summarize(regime_type = first(wth_regime), 
            start_year = first(year), 
            end_year = last(year), 
            mean_input = mean.new(gender_inclusion_input_mean))

table3a <- table3 %>%
  slice_max(order_by = mean_input, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the input side" = mean_input)

stargazer(table3a, summary=FALSE)

table3b <- table3 %>%
  slice_min(order_by = mean_input, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the input side" = mean_input)

stargazer(table3b, summary=FALSE)

## Table C3 ##

table4 <- vdem %>%
  group_by(regime_type_id) %>%
  summarize(regime_type = first(wth_regime), 
            start_year = first(year), 
            end_year = last(year), 
            mean_output = mean.new(gender_inclusion_output_mean))

table4a <- table4 %>%
  slice_max(order_by = mean_output, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the output side" = mean_output)

stargazer(table4a, summary=FALSE)

table4b <- table4 %>%
  slice_min(order_by = mean_output, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the output side" = mean_output)

stargazer(table4b, summary=FALSE)

##  Drop NA on DV side ##
vdem <- vdem %>%
  drop_na(v2xpe_exlgender)

vdem2 <- vdem2 %>%
  drop_na(v2xpe_exlgender)

vdem3 <- vdem3 %>%
  drop_na(v2xpe_exlgender)

#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

## Relevel Factor's ##
vdem$wth_regime <- as.factor(vdem$wth_regime)
vdem<- vdem %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem2$wth_regime <- as.factor(vdem2$wth_regime)
vdem2<- vdem2 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem3$wth_regime <- as.factor(vdem3$wth_regime)
vdem3<- vdem3 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

#### Model Building step by step ####
## compare https://strengejacke.github.io/mixed-models-snippets/random-effects-within-between-effects-model.html 

## Growth Curve model

growth_model <- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party + 
                       (1 + year | country_name),
                     data = vdem3, 
                     REML = TRUE) # REML estimation for parameters
summary(growth_model)
isSingular(growth_model, tol = 1e-05) # FALSE

## Complex REWB ##

complexREWB<- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
                     (1 + year | country_name) + (1 + v2xps_party_dm | country_name),
                   data = vdem3, 
                   REML = TRUE) # REML estimation for parameters
isSingular(complexREWB) # FALSE

complexREWB.1 <- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
                        (1 + year | country_name) + (0 + v2xps_party_dm | country_name),
                      data = vdem3,  # assume independence between random slopes and no coveriance
                      REML = TRUE) # REML estimation for parameters
isSingular(complexREWB.1, tol = 1e-05) # FALSE

## Simple REWB model ##

simpleREWB <- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
                     (1 + year | country_name), 
                   data = vdem3,
                   REML = TRUE) # REML estimation for parameters
isSingular(simpleREWB, tol = 1e-05) # FALSE

## Mundlak model ##

mundlak <- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party + v2xps_party_gm + 
                  (1 + year | country_name), 
                data = vdem3,
                REML = TRUE)
isSingular(mundlak, tol = 1e-05) # FALSE

#### Comparison of Models and check suitable model ####

tab_model(growth_model, complexREWB, complexREWB.1, simpleREWB, mundlak,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Classical RE", "Complex REWB", "Complex REWB (2)", "Simple REWB", "Mundlak"))

## REWB or simple RE-model ? ##
anova(growth_model, mundlak) # significant improvement of the Mundlak-model over the simple 
# RE-model, indicating that it makes sense to model within- and 
# between-subjects effects

## Simple or Complex REWB? ##
anova(simpleREWB, complexREWB)
# Complex REWB Model: significant improvement of the Complex REWB over the the simple REWB model. 

## REWB model or FE model ? ##

fixed <- felm(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm  | country_name,
              data = vdem3)

tab_model(fixed, growth_model, complexREWB,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("FE-model with ID", "Classical RE", "Complex REWB"))

## fixed or Complex REWB? ##

anova(fixed, complexREWB) # no significant improvement of the complex REWB over the simple 
# REWB-model, indicating that it makes no sense to model random slopes and covariance

#### Residual diagnostics for complex REWB model ####

## Package ## 
library(DHARMa)
citation("DHARMa")

## Calculating scaled residuals ##

simulationOutput <- simulateResiduals(fittedModel = complexREWB, n = 500)
simulationOutput <- recalculateResiduals(simulationOutput, group = vdem3$country_name)

## Plotting the scaled redisuals ##

plot(simulationOutput, quantreg = FALSE) # plot leads to the conclusion that there is no big problem with outliers or uniformity

## Goodness-of-fit tests ##

# Uniformity: ok
testUniformity(simulationOutput = simulationOutput)

# Over/underdispersion: no problem 
testDispersion(simulationOutput)

# Heteroscedasticity: no problem
plot(simulationOutput, quantreg = FALSE) 

rm(mundlak, simpleREWB, complexREWB, complexREWB.1, fixed, growth_model, curWarnings, simulationOutput)

############################################################################################################
############################################################################################################

#### Main Models, M1 - M4 ####

m1<- lmer(v2xpe_exlgender ~ year + wth_regime + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)

m3<- lmer(v2xpe_exlgender ~ year + wth_regime + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)


#### Extracting findings and plotting effects ####

tab_model(m1, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4"))

library(texreg)
texreg(list(m1, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "No Party",
                             "One-Party Regime", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C4", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")


m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       wth_regimemonarchy = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party")) %>%
  filter(term!= "Intercept" & term!="Year")


m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       wth_regimemonarchy = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Multi-party" | term=="Monarchy" | term == "No-party" | term =="One-party")


m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.1),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
} 
ggsave("calculations/Output/robustness_wth/Figure3.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m3)


####  Models with Party Institutionalization as Indepedent Variable ####


m5<- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)


m7 <- lmer(v2xpe_exlgender ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)
isSingular(m7, tol = 1e-05)


#### Texreg Findings ####

tab_model(m5, m7,  
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 7"))

texreg(list(m5, m7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "One-party",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C5", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")


m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       wth_regimemonarchy = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeone-party" = "One-party", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       wth_regimemonarchy = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeone-party" = "One-party", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "One-Party" | term=="Monarchy" | term == "Multi-Party" | 
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")


m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/robustness_wth/Figure5.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m7)


############################################################################################################
############################################################################################################

#### Input Side ####

##load data ##

vdem <- readRDS("calculations/data/vdem_wth_merged.rds")

#### Prepare Data for analysis ####
library(scales)

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(1,0))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(1,0))
summary(vdem$gender_inclusion_output_mean)

vdem <- vdem %>%
  drop_na(cown, year, wth_regime, gender_inclusion_input_mean)

vdem2 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)

#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

## Relevel Factor's ##
vdem$wth_regime <- as.factor(vdem$wth_regime)
vdem<- vdem %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem2$wth_regime <- as.factor(vdem2$wth_regime)
vdem2<- vdem2 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem3$wth_regime <- as.factor(vdem3$wth_regime)
vdem3<- vdem3 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

#### Main Models, M1 - M4 ####

m1<- lmer(gender_inclusion_input_mean ~ year + wth_regime + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)


m3<- lmer(gender_inclusion_input_mean ~ year + wth_regime + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)
isSingular(m3, tol = 1e-05)

#### Extracting findings and plotting effects ####

tab_model(m1, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 3"))

library(texreg)
texreg(list(m1, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "No-Party",
                             "One-Party Regime", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C6", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party")) %>%
  filter(term!= "Intercept" & term!="Year")


m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Multi-party" | term=="Monarchy" | term == "No-party" | term =="One-party")



m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.1),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
} 
ggsave("calculations/Output/robustness_wth/Figure3_Input.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m3)

####  Models with Party Institutionalization as Indepedent Variable ####

## Relevel Factors ##
vdem3$wth_regime <- as.factor(vdem3$wth_regime)
vdem3 <- vdem3 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))


m5<- lmer(gender_inclusion_input_mean ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)

m7 <- lmer(gender_inclusion_input_mean ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)
isSingular(m7, tol = 1e-05)


#### Texreg Findings ####

tab_model(m5,  m7,  
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 7"))

texreg(list(m5, m7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "One-Party Regime",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C7", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")


m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")


m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term=="Monarchy" | term == "Multi-party" | term =="One-party" |
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")


m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/robustness_wth/Figure5_Input.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m7)

##############################################################################################
################################################################################################
################################################################################################

#### Output Side ####

##load data ##

vdem <- readRDS("calculations/data/vdem_wth_merged.rds")

#### Prepare Data for analysis ####
library(scales)

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(1,0))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(1,0))
summary(vdem$gender_inclusion_output_mean)

vdem <- vdem %>%
  drop_na(cown, year, wth_regime, gender_inclusion_output_mean)

vdem2 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, wth_regime, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)

#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

## Relevel Factor's ##
vdem$wth_regime <- as.factor(vdem$wth_regime)
vdem<- vdem %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem2$wth_regime <- as.factor(vdem2$wth_regime)
vdem2<- vdem2 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

## Relevel Factor's ##
vdem3$wth_regime <- as.factor(vdem3$wth_regime)
vdem3<- vdem3 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))

#### Main Models, M1 - M4 ####

m1<- lmer(gender_inclusion_output_mean ~ year + wth_regime + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)


m3<- lmer(gender_inclusion_output_mean ~ year + wth_regime + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)


#### Extracting findings and plotting effects ####

tab_model(m1, m3,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 3"))

library(texreg)
texreg(list(m1, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "No-Party Regime",
                             "One-Party Regime", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C8", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")


m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party")) %>%
  filter(term!= "Intercept" & term!="Year")

m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Multi-party" | term=="Monarchy" | term == "One-party" | term =="No-party")


m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.1),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
} 
ggsave("calculations/Output/robustness_wth/Figure3_Output.pdf", height = 14, width = 24, units= c("cm"))


## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m3)


####  Models with Party Institutionalization as Indepedent Variable ####

## Relevel Factors ##
vdem3$wth_regime <- as.factor(vdem3$wth_regime)
vdem3 <- vdem3 %>% 
  mutate(wth_regime = relevel(wth_regime, "military"))


m5<- lmer(gender_inclusion_output_mean ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)

m7 <- lmer(gender_inclusion_output_mean ~ year + wth_regime + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)
isSingular(m7, tol = 1e-05)


#### Texreg Findings ####

tab_model(m5, m7,  
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 7"))

texreg(list(m5, m7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Monarchy",
                             "Multi-Party Regime",
                             "One-Party Regime",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C9", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")


m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "wth_regimemonarchy" = "Monarchy", 
                       "wth_regimemulti-party" = "Multi-party", 
                       "wth_regimeno-party" = "No-party", 
                       "wth_regimeone-party" = "One-party",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Multi-party" | term=="Monarchy" | term == "One-party" | term =="No-party" |
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")

m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/robustness_wth/Figure5_Output.pdf", height = 14, width = 24, units= c("cm"))


## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m7)
